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A NEW ITERATION FOR LOCATING EQUILIBRIUM POINTS 

IN NONLINEAR SYSTEMS 
Clarence Cantor, Systems Division 

and 

Fawzi P. Emad, Department of Electrical Engineering 

Abstract 

A discrete nth order nonlinear dynamic system, or an iteration for a nonlinear set 
of n algebraic equations in n unknowns, can be represented by 

x k.i = f ( x k) 

where x k is an nth order state vector. The problem is to locate the equilibrium points 
of the system, give some approximation of these points, when the original iteration is 
divergent or only slowly convergent. Newton's method solves the problem in general, 
but requires f(x) to be available in analytic form and also requires extensive compu- 
tations of partial derivatives. Steffensen's iteration, utilizing sets of n + 2 iterates 
obtained from the original iteration, can also solve the problem while avoiding the dif- 
ficulties in Newton's method. Steffensen's iteration is based on a linearization of the 
system about an equilibrium point x e , in the form 

x k+i ■ x « = A ( x k' x .) • 

However, Steffensen's iteration breaks down when the A matrix corresponding to a set 
of iterates has an eigenvalue equal to one, or when the matrix cannot be determined. 
The proposed iteration is equivalent to Steffensen's iteration in the general non- 
singular case, but it can be readily extended in a meaningful way in a large class of 
singular cases where Steffensen's iteration breaks down. Also, the. proposed method 
should produce better results in near singular cases where Steffensen's iteration can 

cause discrepancies due to numerical errors. 

iii I 


fMCEPINO PAGE BLANK NOT..EILM££U 

i 5 


CONTENTS 

Page 

ABSTRACT . . iii 

INTRODUCTION . . . 1 

PROPOSED NEW METHOD 6 

EXAMPLE * 21 

CONCLUSIONS 23 

REFERENCES . „ 24 


v 


INTRODUCTION 


A large class of problems in nonlinear systems involves determining the solutions 
(equilibrium points) of the vector equation, 

x - f(x) . (1) 

For example, in a discrete nonlinear dynamic system represented by 

x k + i = f ( x k) (2) 

where x k is an nth order state vector, it is often necessary if not mandatory to determine 
the equilibrium points x e (there may be more than one), where x g is defined by 

x e = f ( x e) • (3) 

Equation 2 can also represent any iteration scheme involving n unknowns to be 
determined „ 

If an approximation x Q of an equilibrium point is known, and if the process defined 
by (2) is convergent (asymptotically stable) in a region containing x Q , then the equilibrium 
point x e can be obtained by repeated application of (2). However, this is not possible if 
the process is divergent, or impractical if the process is only slowly convergent. 

Newton* s method [1] , [2] is a powerful method that can be used in general to solve 
equation 1 (or 2), even when the iteration defined by equation 2 is divergent. Let 

g(x) = f(x) - x v (4) 

Then the equation to be solved is 

g(x) = 0 : . ( 5 ) 
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If x Q is an approximate solution of (5)' (or (1)), and x e is the exact solution, then (5) 
can be linearized about x = x 0 to yield 

g( X o) + j( X o)( X e" X o) = 0 (6) 

where J (x Q ) is the Jacobian of g( x) evaluated at x = x Q . 

Assuming J -1 (x Q ) exists, we can solve (6) for an approximation of x e , 

x e * x o " J -1 ( x o)g( x o) * (7) 

Since this is only an approximate solution of x e (although hopefully a better approxi- 
mation than x Q ) , we can treat (7) as an iteration, namely 

X 1 = x » “ J" 1 ( X oM X o) and x k+l = x k ” J _1 ( x k)«( X k) < 8 > 

which is Newton’s iteration for a system of equations. 

In terms of (1) , we have 

s( x k) = f ( x k)‘ x k and J( x k) = A ( x k) - 1 

where A(x k ) is the Jacobian of f(x) at x = x k . Hence, (8) can be written as 

x k*i = x k - ( A ( x k)- 0 1 ( f ( x J- x k) • (9) 

Equation (8) (or (9)) will converge to x e if x Q is sufficiently close to x e . One of the 
difficulties in Newton’s method is the need to calculate J (x k ), with its n 2 partial deriva- 
tives, at each x k , which can require rathpr extensive calculations. This is alleviated 
in the modified Newton’s method which in general is not as rapidly convergent. Here 
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J (x Q 'j is retained throughout the iteration, namely 


x k+i " x k “ J _1 ( x o) g( x k) * (10) 

Another limitation of Newton's method is that g(x) (or f(x)) must be available in 
analytic form. Thus a computer simulation of a discrete control system x k+1 = f (x k ] , 
where f(x) is not available in analytic form, and only sampled output x k is available, 
cannot be treated directly by Newton's method. 

Steffensen's iteration [3] is another method for solving (2) which is similar to 
Newton's method. It has the advantage of not requiring the calculation of partial deriva- 
tives. Also, it does not require that f(x) be available in analytic form. It utilizes sets 
of n + 2 iterates obtained either from the system equation (2), or from a computer 
simulation of the system. 

Assume that (2) is approximately linear in a region surrounding the equilibrium 
point x e . This linearity can be expressed as 


X. ,, - X 

k+1 e 


- A K^ X e) 


( 11 ) 


where A represents the Jacobian of f(x) at x e . 
It is easy to show then that 


Ax,.,. = x, - x, - A(x k+1 ~x k ) 


k+l “ k+2 k+1 


or 


Ax k+1 - AAx k , 


( 12 ) 
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We define nx n matrices X 0 , AX^, 'and AX x as follows. 

X 0 = ( X 0 X 1 X 2 

= (Ax 0 A X 1 •■•-Ax„_ 1 ) 

AXj = (Axj Ax . 2 * • * Ax n ) 

Then (12) yields 

AX x = AAX 0 or A = AX t (AX 0 ) _1 
assuming (aX q ) _1 exists. 

Utilizing (11), withk = 0, we can solve for x e as follows. 

(A-I)x e = Ax 0 - x 0 - (x t -x Q ) = (A-I)x 0 -Ax 0 

or 

X e = X 0 " ( A ” 1 )" 1 Ax 0 • 

Using the expression for A of (14), we obtain after some manipulation, 

(A-I r 1 = AX 0 (AXj-AX,,)- 1 

and 

X. = *0 - • 


(13) 


(14) 


(15) 


(16) 


Since the system is not really linear, (16) yields only an approximation of x e , but 
one which hopefully is closer to x e than the starting point x Q . Calling this new approxi- 
mation x^ 1 ) , we have 

x 0 (1) = (17) 
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We can repeat this process using x Q (1) as the Starting point for a new set of n + 2 
iterates, and apply (17) again to obtain and so forth. This yields the sequence 

x Q (0) , x Q (1) , x 0 ^ 2) , ° • • o This process is Steffensen’s iteration. Hopefully, the se- 
quence x Q (0) , x Q (1) , x 0 < 2 >, • • will converge to x e . In general, this convergence will 

occur if the first starting point x Q (0) is close enough to x g , even when the original itera- 
tion of equation 2 diverges. 

One of the difficulties in Steffensen's iteration occurs when the matrix (AX x - AX Q ) 
is singular or nearly singular. When it is singular, there is nothing one can do to con- 
tinue the process except perhaps introduce an arbitrary perturbation in the starting point 
x 0 (k) and hope that the new set of n + 2 iterates yields a non-singular (AX x - AX 0 ). When 
the matrix (AX 1 - AX Q ) is nearly singular, numerical errors can produce a large dis- 
crepancy in the resulting new approximation of x e . 

The singularity of (AX x - AX 0 ) can arise in two ways. Since AX x = AAX Q , we have 

- AX„ - (A -I) AX 0 . (18) 

Thus AX : - AX Q will be singular if either (A- 1) or AX Q is singular. 

The proposed method to be described next provides a systematic and meaningful 
extension of the iteration when (AX x - AX Q ) is singular or nearly singular. It avoids the 
numerical errors associated with inverting a near singular matrix. In the general 
case, when (AX x - AX Q ) is non-singular, the proposed algorithm yields the same results 
as Steffensen’s iteration with an equivalent amount of computations. It has one advan- 
tage here also in that after an iteration, the A matrix can be obtained with less additional 
calculations than would be the case in Steffensen’s iteration. 
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PROPOSED NEW METHOD 

First we will consider the general case when (AX 1 - AX Q ) is non-singular. This 
implies that neither (A- 1) nor AX 0 is singular. 

Consider the matrix 


AX' = ( Ax n Ax •••••■Ax Ax ) 

0 .v 0 1 n-1 n / 


This matrix is obviously singular since it contains n + 1 columns and only n rows. 
Hence, the last column, Ax n , can be written as a linear combination of the first n 
columns, namely 


Now 


Ax c n Ax_ + c, Ax, +••■•+ c , Ax , 

n 0 0 1 1 n-1 n-1 


Ax : £ x... - x 


i - ( x i + i ~ x .) - ( x i - x .) 


- (A-I> ( x A - x e ) 


utilizing (11). Thus we can write (20) as 


(A-I)( x „- x «) - (A-!) [ c o( x o- x .) + c 1 ( x i- x ,) + "" +C „-i ( x n-! “ X e)] • <22 

Premultiplying by (A -I)" 1 (which is assumed to exist in this general case) yields 


V " x . - c o( x o- x .) + c i (*i _x «) + :‘ + c „-, ( X „-1 - x .) • 
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Solving for x e yields 


Assuming 


c x + c, x, + • • • + c , x - x 

00 11 n~l n-1 n 


n~ 1 

L 

i—0 


(24) 


C i " 1 


n-1 

ZN * 1 


j=0 


(which we will prove), (24) represents our next estimate of x e which we will call x^ 1 ) . 
Thus 



x. + 


0 0 


c l x l 


+ 


+ c 


n-1 n-1 


” X 



<2 


&) 


To solve for the c., we go back to (20) which can be written as 


- (Ax 0 Ax 1 ••• Ax n _i) 


The n x n matrix is simply AX Q whose inverse is assumed to exist. Hence 


(AX 0 )" 1 Ax n - ( AX Q )“ 1 (x n+1 -x n ) . (27) 




7 



The combination of (25) and (27) represents the new algorithm for determining 
x Q (1) , x q (2) , etc. Equation (27) is used first to determine the constants c i which are 
then used in (25) to obtain x 0 (1) . Then x 0 (1) is used as the starting point to obtain a new 
set of n + 2 iterates from (2), and the process is repeated. This algorithm for the ap- 
proximation of x e has been derived using the same equation (21) that yields Steffensen’s 
iteration, and thus gives identical results in the general non-singular case. The amount 
of calculations in each algorithm is about the same. The new algorithm has the advantage 
of yielding the A matrix, when that is desired, with fewer additional calculations than 
Steffensen’s iteration. Since ( AX 0 ) _1 is calculated in the new algorithm, A = (AX x ) AX 0 _1 
is determined from one additional matrix multiplication. In Steffensen’s iteration 
(equation 17), we determine the matrix AX 0 (Z\X 1 - AX 0 ) _1 = (A- 1) -1 , so that a matrix 
inversion is required before extracting A. The greatest advantage of the new algorithm, 
however, is that the same general form is used to extend the iteration in singular or 
near singular cases, where Steffensen’s iteration breaks down or causes large errors. 
This will be discussed later. 

We still have to prove that 

n— 1 

£> / 1 • 

i=0 

in order for the new algorithm to be valid. 

Proof: Assume 


n-1 



i=0 
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Then (20) can be written as 


n— 1 

= C 0 AX 0 +C 1 AX 1 + + C n-l AX n-l 


( 28 ) 


i=0 


Using (12), we can write Ax n as 


Ax - A Ax , - A 2 Ax „ - - A n Ax. 

n n~l n~2 0 


(29) 


Utilizing (29) in (28) , we obtain 

c Q (A n -l)Ax 0 + c x (A n_1 - i) Ax 1 +•■•• + + c n-1 (A- 1) Ax n-1 = 0 . (30) 

Each term in parenthesis in (30) contains the factor (A -I) whose inverse exists 
in the general case. Then multiplying equation 30 by (A - 1) - 1 yields 

c Q (A n_1 +A n " 2 + • • ••+ A+ i) Ax 0 + c 1 ( A n “ 2 + A n " 3 + .** - A + I ) Ax x 

+ + c n _ 2 (A + I) x n _ 2 + c n _ x Ax n-1 = 0 . ( 3i) 

Performing the indicated matrix multiplications, using (12), yields 

C 0 Ax 0 + ( C 0 + C l) Ax i + ( c 0 + c i + + c n-l) Ax n-1 = 0 

or 

Ax „-1 = - c 0 Ax 0 - ( c 0 +C l) Ax i - ( C 0 +C l + •"' +C n- 2 ) ■ (32) 

Equation (32) implies that the last column of AX 0 is a linear combination of the 
first n - 1 columns. This implies that det AX Q = 0 which contradicts the fact that 
( AX q ) “ 1 exists in the general case. Q.E.D. 
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Thus 


n~l 

Z> *■ 1 


i=0 


and the algorithm of (25) and (27) is valid. 

We will now extend the algorithm to cover the case of | AX 0 = 0 and | A— 1| f 0. 
In all of the subsequent arguments, we assume that Ax Q ^ 0, because if it were the 
problem is solved trivially, i.e. A x Q = 0 => Ax 0 = Ax x = 0 which implies that 

x Q = x 1 = x 2 = = x n+1 and we are at an equilibrium point already. 

Let r be the rank of the matrix (AX q ). Then the first r columns of (AX q ) are 
linearly independent. 

Proof: Let k be the maximum number of consecutive columns, starting with the 
first, that are linearly independent. Then 1 A k < r < n. Then the k + 1 column must 
be expressable as a linear combination of the first k columns, or 


Ax, - a 
k 


0 Ax 0 + a i Ax ! + • ••*•+ a k-i Ax k-1 


(33) 


We can then express Ax k+1 as 


Ax k+1 = AAx k - a 0 Ax i + a l Ax 2 + + a k-l Ax k 


^k+l = a 0 AX 1 + a i ^2 + — " + a k~2 ^k-l + a k-l ( a 0 Ax 0 + a l **1 + + a k-l ^k-l) 


or 


Ax k + 1 " b o Ax o + b i Ax i + *■ b k-l ^k-l • 


(34) 
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Thus the k + 2 column is also expressable as a linear combination of the first k 


We can express the k + 3 column, Ax. 2 , as 


Ax k+2 “ A ^kH " b 0 Ax i 4 b l AX 2 4 + b k-2 ^k-l 


4 b k-i ( a o Ax o 4 a i Ax i 4 4 a k-i Ax k-l) 


Thus the k + 3 column is again a linear combination of the first k columns. Con- 
tinuing the process will show that every column after the first k columns is a linear 
combination of the first k columns. Hence the rank of AX Q must equal k or k - r. 

Since the first r columns are linearly independent, we can express the r + 1 
column as a linear combination of the first r columns. Thus 


Ax - c. Ax. + c, Ax. + + c , Ax . 

r 0 0 11 r-1 r-1 


From (21) , we have 


Ax i - (A - I) (x. - x e ) . 


Then (36) can be written as 


(A-I) (x r -x e ) - (A- I) [c 0 (x 0 - x e ) + * + c r ^ (x r _ x -xj] 


Premiiltiplying by (A-I) 1 and solving for x e ,we obtain 


C 0 X 0 + C 1 X 1 4 4c r-l X r-1 ~ X r 
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Since this represents our next estimate of x e , we denote it by x n (1) , or 


x 

x o 


C 0 X 0 + C 1 X 1 + + C r-l X r-l- X r 

r-1 


Z> - 

V i= 0 / 


( 38 ) 


Note that this is exactly the same form as (25) for the general non-singular case, 
except that n is replaced by r . We can prove, that 

, i 

i-0 


in the same way that we proved that 


n-j 


/ l 


i=0 


in the non- singular case. Simply replace n by r in the latter proof. 

To complete the algorithm of (38), we must solve for c., i = 0, 1, • • • r - 1. If 
we take the dot product of (36) with Ax Q , Ax , • « • , and Ax r _ 1 , we obtain 


c o 


Ax o Ax r ~ ” 
Ax, • Ax 

1 


1 r 

^r-JL 


Ax , ••Ax 

r-1 r 


(39) 
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where 


Ax o "^o 
Ax, • • Ax„ 


Axq • Ax x 
Axj * Z\x 1 


AX 0 * AX r - 1 

Ax 1 • Ax f _ 1 


Ax 


r-l 


Ax, 


Ax 


Ax, 


Ax 


r-l 


Ax 


r-l 


is the Gramian corresponding to the first r columns of AX . Since these columns are 


linearly independent, G r f 0 and G r 1 exists. Then, 



Ax 0 ••Ax r 
Ax x • Ax r 



(40) 


Equations (38) and (40) constitute the algorithm in the case where A X Q is singular 
with rank r and (A- 1) is non-singular. Note that these equations can also serve as the 
algorithm for the general non-singular case by letting r = n. Then (38) becomes the 
same as (25) and (40) reduces to (27) since 


and 


G 

n 


= (®o) T ®o 


Ax n • Ax 

0 n 

Z^c.. • Ax 

1 n 




= (^») T& n • 



Then 



o; 1 K)\ = ^‘KT 4 * o^n 


X” 1 Ax 

0 n 


( 43 ) 


which is the same as (27). 

We will now show that the algorithm of (38) and (40) can be used meaningfully in 
many cases when (A -I) is singular. The only restriction on its use is that the null 
spaces of (A -I) and (A -I) 2 be equal. This is equivalent to the condition that (A -I) is 
fully degenerate i.e., the number of independent eigenvectors V i such that (A- 1) V i = 0 
equals the multiplicity of the zero eigenvalue of (A -I). The proof of this equivalence 
has been omitted for the sake of brevity. 

Again we assume that Ax 0 f 0 for if it were, we would already be at an equilibrium 
point. We have 


- (Ax 0 Ax, Ax n-1 ) 

= (A-I) [(x 6 -xJ( Xl -x e ) (x^j-xjl . (44) 

Then 

| A - 1 1 = 0 => |AX 0 | = 0 : . 

Let Vv, V 2 , • • • V k be k independent eigenvectors of (A- 1) associated with the zero 
eigenvalue of (A -I), where k is the maximum number of independent eigenvectors. Then 


(A-iyy. = 0 i = 1, 2 , k . 


(45) 
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These k eigenvectors thus span the null space of (A - 1) . 


Let r equal the rank of AX Q where r < n. Then the first r columns of AX Q are 
linearly independent as proved previously. The r + 1 column can then be expressed 
as a linear combination of the first r columns, or 


Ax - c. Ax rt + c, Ax, + • • • + c , Ax , 

r 0 0 11 r-1 r-1 

Using equation 21, we obtain 


(46) 


(Arl) (x r -x e ) - (A-I) [c 0 (x 0 -x e )+c 1 (x 1 -x e )+ + V 1 (x r . 1 -x e )] . (47) 

Let V be some vector, as yet undefined, in the null space of (A-I). Then V can be 
written as a linear combination of the k eigenvectors V i5 or 

V = b i V i + b a V 2 + .-. +b k V k . (48) 

Then 

(A-I) V = 0 . (49) 

Using (49) in (47), we obtain 

(A-!) [ c o( x o“ x ,“ V ) +c i( x i“ x e“ V ) 

+ + C r 1 ( x r 1 ■ X „ ■ V ) ■ ( X r “ X * ■ V )] = 0 • (50) 

Each vector (x. - x e ) can be divided into two components, one orthogonal to the 
null space of (A-I) and the other in the null space of (A-I). Denoting these components 
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by x.' a and x.' b respectively, we have 


& i i / 

X. - X = X. + X., 
i e la lb 


i - 0 , 1 , • - • r 


(51) 


By definition, 


(A ~ I) x ' = 0 . 


(52) 


Substituting equation 51 in equation 50, we obtain 


(A-I)jcx' + c x ' + +c x ' - x ' +c(xl-v) 

v ' L 0 0a 1 la r-1 r-1, a ra 0 \ Ob / 


+ + C r-1 ( X t-l,b ^ V ) " ( X r'b - V )] = 0 • 


(53) 


We can select V by proper choice of the constants b i of equation 48, so that 


V = 


C X ' + C , X ' +••••■+ C ■ X ' . - X ' 

0 Ob 1 lb r-1 r-1, b rb 

r-1 


L 

i=0 


(54) 


It can be shown that 


' r-1 \ 

h Ci ) * 1 ’ 


i-0 


provided that the null spaces of (A- 1) and (A-I) 2 are equal. The proof has been omitted 


here for the sake of brevity . 
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Using (54) in (53) we obtain 


(A - I) [c. x ' + c, x,' + • ••• •+ c x ' - x ' I - 0 

v J L 0 Oa 1 la r-1 r-1, a ra J 


(55) 


Let 


= c n X ' + c. x' + 

0 Oa 1 la 


+ C , X , 
r“l r-l,i 


“ X 


r a 


(56) 


Then (A -I) z = 0. 

We can easily prove that z = 0. 

Proof: Assume z f 0. Then z is orthogonal to the null space of (A -I) since each 
of its components is orthogonal to this null space. This implies (A- 1) z f 0 which is 
a contradiction. Thus 


z - c. x ' + c, x/ + c . x ' - x ' - O' 

0 Oa 1 la r-1 r-l,a ra 


(57) 


Equation 54 implies that 


s( x o'b- v ) + c i( x i; - v ) + "- ■ +c rt(*A.b- v ) - ( x ,b- v ) = o- 


(58) 


Adding (57) and (58), and using the identity of (51), we obtain 


C 0 ( X 0 " X e " V ) + C 1 ( X 1 - X e " V ) + ’ ‘ - + C r-1 ( X r-1 " X e “ V ) " ( X t” X e" V ) ' 0 V (59) 

We note that the point ( x e + V ) is an equilibrium point, which we will denote by x e ' , 
since 


A(x;-x e ) = AV = V = (x; -xj . 


(60) 
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We can solve for x e ' = x e + V from equation 59 to obtain 

Vl X r-l" X r 

(61) 

- 1 

Equation 61 indicates that we can obtain an approximation of an equilibrium point 
x e ' , which in general is not the same as x e , when (A -I) is singular. Of course this is 
the best that any iteration scheme can do. If the system were truly linear with | A - 1] - - 0, 
then the next set of iterates starting at x e ' would remain stationary at x fe ' . The ques- 
tion arises, in a nonlinear system with an isolated equilibrium point, can the A matrix 
corresponding to a set of iterates be such that |A-l| = 0? We know that this can be 
true at the isolated equilibrium point itself, but of course this is a trivial case since 
this would mean we were already at the equilibrium point. Without answering the 
question (which may be impossible to answer), we can at least say that the matrix (A -I) 
can become near singular, as we get closer to an equilibrium point whose (A- 1) matrix 
is singular. In practical computations the difference between a singular matrix and 
near-singular matrix can be academic since the numerical errors associated with 
inverting a near singular matrix can cause large discrepancies. It is likely that treat- 
ing the matrix as singular, when - e < | A - 1 1 < e , will yield more convergent results in 
many cases. Thus equation 61 can serve as the algorithm for the next estimate of the 
equilibrium point, when (A- 1) is singular or near singular, or 




This is identical with (38), the algorithm for obtaining an estimate of x e when 
|AX 0 | = 0 and j A — I J f 0. Thus without knowing whether or not jA-ll = 0 we can use 
the same algorithm for determining x Q (1) whenAX 0 is singular or near singular. The 
c. in (62) are derived from the same equation that yielded the c L for (38). Hence 


c o 

C 1 


Ax q - Ax,. 
Ax 1 • Ax r 

• • 

II 

O 

i 

»-» 

•# ' 

* • 
• ■ 

r 

Hr 
*1 • • 

S’ 

L: 


(63) 


We have thus established that the algorithm of (38) and (40) (or (62) and (63)) is 
valid for the case of | AX 0 1 = 0 regardless of whether or not | A— 1| = 0, provided that 
if 1 A — 1| = 0 the null spaces of (A-I) and(A-I) 2 are equal. We will show that the 
algorithm can also be used effectively in the case of near singularity of AX 0 in the sense 
that the c. obtained from (40) (or (63)) will minimize the norm of y, where 


y = c Ax + c, Ax, + 


■ + c 


r-l 


Ax^j-AXr 


(64) 


If AX 0 were truly singular and of rank r , then the selection of c i as per equation 40 
(or 63) would make y identically zero . However assume that AX 0 is nearly singular 
and that it is desired to treat AX 0 as though it were of rank r where r is chosen to be 
the smallest integer satisfying a predetermined condition, 


- e < 



< e 


e > 0 : . 


In other words, r is chosen as the smallest integer such that the r + 1 column can 
"almost’ 1 be expressed as a linear combination of the first r columns. We would like 
the difference between the selected linear combination and the r + 1 column, Ax r , to 
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be as "small" as possible. This difference vector is y. We will consider the choice 
of c. to be optimum if ||y|| is minimized where 



We will now derive a formula for the c. that minimizes ||y || (or ||y|| 2 ). 

II y II 2 = yy = ( c 0 Ax o + c i Ax i + + c r -i Ax r -i _Ax r ) 

" ( C 0 AX 0 +C 1 AX 1 + •"•' +C r-1 AX r-l ~ Ax r) * 

Setting 

^ lly I ! 2 _ n . 

dc. - 0 > 

X 

i = 0, 1 , • • • , r - 1 , we obtain r equations of the form, 

0 ; - 2(c n Ax. + c Ax, c Ax -Ax \ ••Ax- . 

V O 0 11 r— 1 r-1 r / x 

This yields the matrix equation 


Ax q * Ax 0 Ax q • Ax 1 Ax q • Ax r-1 


c 0 


Ax 0 * Ax r 

Ax. * Ax. Ax, "Ax, Ax, *-Ax , 

10 1 1 1 r-l 


c i 

• > 

_ 

Ax, *-Ax 

Ax r _ 1 • Ax Q Ax r-1 •■Ax 1 Ax r-1 ••Ax r _ 1 


•: * 

l r ~L 


Ax • Ax 

r-l r 


( 65 ) 


( 66 ) 


(67) 
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or 


c o 


Ax 0 • Ax r 

C 1 


Ax x • Ax r 

• - 

= G -1 


• 

l r ~L 


Ax • • Ax 

r-l r 


( 68 ) 


This is identical to (40) (or (63)) so that selection of the c L in accordance with the 


algorithm will minimize y and thus yield a "best" fit in the case where AX 0 is nearly 


singular. In the above derivation, it was not proven that ||y is a minimum rather than 


a maximum. However the fact that ||y{] - 0 as 


r+1 


0 indicates that i|y|| must be a 


minimum for the c. chosen as per (68), 


EXAMPLE 


To illustrate some of the points that have been discussed, a simple example will 
be presented in which the initial starting point x Q results in singular matrices AX 0 and 
( AXj - AX q ) . Obviously, Steffensen's iteration would break down in this instance. Let 

x =' y + x 2 y and y - x + xy 2 . (69) 

This obviously has an equilibrium point at (0, 0). However, we will pretend that 
we only have some approximation which happens to be (0.1, 0.1). Forming n + 2 
iterates from (69), we obtain 


- m 


r o : . ioii 


r. 1020301 


r. 1030921 


x ° Lo. ij 

x i 

LO. 101 J 

% 2 

L. 102030 j 

X 3 - 

L. 103092 J 

(70) 

ii 

o 

<5 

rovooiri 

H = 

r ooio3oi 

Ax 2 ■ - 

r 0010621 



LO.OOlJ 

1 001030 J 

L 00106 2 J 


(71) 
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We can see that AX 0 and (AX x - AX 0 ) are singular. Using the algorithm of (38) 
and (40), we obtain 

. 0666671 

. 066667 J ’ ( 72 ) 

Using x^ 1 ) as the starting point for a new set of iterates, we obtain 


0666671 
|_. 066667 J 

x i “ 

r. 0669631 
[_. 066963J 

X 2 - 

f. 0672631 
L. 067 263 J 

(73) 

Ax 0 = 

r. 0002961 
L 000295 J 

Ax, = 

0003001 
000300 J 

• ' 

(74) 


It is not necessary to go beyond x 2 since we recognize that AX 0 is singular. Using 
(38) and (40) once more, we obtain 

f. 0447571 

c - 1-01351 and x - * (75) 

o o l_. 0447 57 J {/t>) 

Repeating will show that the sequence x Q (1 ) } x 0 (2) , « • • converges to (0, 0). The 
relatively slow convergence is due to the fact that in our example, A has an eigenvalue 
that approaches one as we approach the equilibrium point. Any initial condition of the 
form ( a, a) is along the eigenvector corresponding to this eigenvector, resulting in a 
singular AX Q . 

This example also has an eigenvalue that is approximately equal to minus one near 
the equilibrium point. If we select an initial condition of the form (a, -a) , which is 
along the eigenvector corresponding to this eigenvalue, AX Q will still be singular but 
the convergence of the algorithm will be much faster. For example, starting at 


c Q - 1.03 and x Q (1) - 
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(0.1, -0.1) yields 


' o . r 

x i 

1011 

X 0 

‘ .102030" 

_-0. 1. 

1 

i 

O 

L— 

2 

L-. 102030. 


.2030301 
203030 J 

The algorithm of (38) and (40) yields 


Ax, 


20 T 
.201 


( 76 ) 


(77) 


c Q - -1.01001 and x Q (1) 


" .0000005 
.-.0000005 


(78) 


which already is extremely close to the equilibrium point. 

CONCLUSIONS 

The algorithm of (38) and (40) is a general method for obtaining the equilibrium 
points of the vector equation x k+1 = f (x k ). When the matrix (AX x - AX 0 ) is non- 
singular, the algorithm yields results identical to Steffensen's iteration with about the 
same amount of computations. The algorithm has one advantage in this case in that it 
enables the determination of the A matrix, when that is desired, with fewer additional 
calculations than does Steffensen's iteration. When the matrix (AX t - AX 0 ) is singular, 
the algorithm yields meaningful results whereas Steffensen's iteration breaks down in 
this case. Finally when the matrix (AX -AX Q ) is nearly singular, (or equivalently when 
AX Q is nearly singular), the algorithm yields good results by treating AX Q as singular. 
This avoids the numerical errors in inverting a near singular matrix. Thus the al- 
gorithm developed in this paper is a more comprehensive method than Steffensen's 
iteration and includes the latter as a special case. 
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